################################################################
#####                 Replication File                       ###
##### The Pacifying Effects of Local Religious Institutions  ###
#####    Alexander De Juan, Jan Pierskalla, Johannes Vüllers ###
################################################################

setwd("~/Dropbox/GIGA/ReligionViolence/")


rm(list=ls(all=TRUE))

library(foreign)
library(lme4)
library(scales)
library(ggplot2)
library(memisc)
library(Zelig)
library(MASS)
library(pastecs)
library(xtable)

# SET PATHS AND LOAD DATA
data <- read.dta("rel_inst.dta")


#Substantive Effect
baseline <- glm(MassFight ~ logdensvil+povrateksvil+natres+urban+golkar1+pdip1+ covyredvil+ wgcovegvil+ wgcovrgvil+ ethfractvil+ relfractvil+ lpop+ ldist+ unemployrate+ DistPolice+ ltotal+ rpol+ relf, family=binomial(link="logit"), data=data)
summary(baseline)

means <- coef(baseline)
var <- vcov(baseline)
draws <- mvrnorm(n=1000, means, var)
concvec <- seq(min(baseline$data$ltotal,na.rm=T), max(baseline$data$ltotal,na.rm=T), by=((max(baseline$data$ltotal,na.rm=T)-min(baseline$data$ltotal,na.rm=T))/(500-1)))

meanvec <- lowervec <- highervec <- NULL
for(i in 1:length(concvec)){
  out <- draws[,1] + mean(data$logdensvil,na.rm=T)*draws[,2]+mean(data$povrateksvil,na.rm=T)*draws[,3]+0*draws[,4]+0*draws[,5]+
    0*draws[,6]+0*draws[,7]+mean(data$covyredvil,na.rm=T)*draws[,8]+
    mean(data$wgcovegvil,na.rm=T)*draws[,9]+ mean(data$wgcovrgvil,na.rm=T)*draws[,10]+concvec[i]*draws[,17]+
    mean(data$ethfractvil,na.rm=T)*draws[,11]+mean(data$relfractvil,na.rm=T)*draws[,12]+
    mean(data$lpop,na.rm=T)*draws[,13]+mean(data$ldist,na.rm=T)*draws[,14]+
    mean(data$unemployrate,na.rm=T)*draws[,15]+mean(data$DistPolice,na.rm=T)*draws[,16]+
    mean(data$rpol,na.rm=T)*draws[,18]+mean(data$relf,na.rm=T)*draws[,19]
  prob <- 1/(1+exp(-out))
  meanhelp <- means[1] + mean(data$logdensvil,na.rm=T)*means[2]+mean(data$povrateksvil,na.rm=T)*means[3]+0*means[4]+0*means[5]+
    0*means[6]+0*means[7]+mean(data$covyredvil,na.rm=T)*means[8]+
    mean(data$wgcovegvil,na.rm=T)*means[9]+ mean(data$wgcovrgvil,na.rm=T)*means[10]+concvec[i]*means[17]+
    mean(data$ethfractvil,na.rm=T)*means[11]+mean(data$relfractvil,na.rm=T)*means[12]+
    mean(data$lpop,na.rm=T)*means[13]+mean(data$ldist,na.rm=T)*means[14]+
    mean(data$unemployrate,na.rm=T)*means[15]+mean(data$DistPolice,na.rm=T)*means[16]+
    mean(data$rpol,na.rm=T)*means[18]+mean(data$relf,na.rm=T)*means[19]
  mean <- 1/(1+exp(-meanhelp))
  lower <- quantile(prob,0.05)
  higher <- quantile(prob,0.95)
  meanvec <- c(meanvec,mean)
  lowervec <- c(lowervec,lower)
  highervec <- c(highervec,higher)
}

pdf("InstitutionEffect.pdf",height=7, width=7*1.62)
plot(concvec,meanvec,type="n",ylim=c(0,0.4), xlab="Religious Institutional Density", ylab="Probability of Mass Fighting")
polygon(c(concvec, rev(concvec)), c(lowervec,rev(highervec)), col="grey80", border=NA)
points(concvec,meanvec,type="l",lty=1)
rug(data$ltotal)
dev.off()
